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Abstract 

In this article we review recent results on the breakup of cylindrical jets of a Newtonian fluid. 
Capillary forces provide the main driving mechanism and our interest is in the description of the 
flow as the jet pinches to form drops. The approach is to describe such topological singularities by 
constructing local (in time and space) similarity solutions from the governing equations. This is 
described for breakup according to the Euler, Stokes or Navier-Stokes equations. It is found that 
slender jet theories can be applied when viscosity is present, but for inviscid jets the local shape of 
the jet at breakup is most likely of a non-slender geometry. Systems of one- dimensional models of 
the governing equations are solved numerically in order to illustrate these differences. 
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1 Introduction 


A perfectly cylindrical column of fluid at rest is an exact solution of the governing equations whether 
viscosity is included or not. When surface tension acts, this configuration is unstable to long waves 
due to capillary instability; any small perturbation with wavelength larger than the undisturbed 
jet radius is unstable. These linear solutions were first given by Rayleigh’s pioneering studies in 
capillary instability, see for example Rayleigh [17], [23], For inviscid jets, for instance, Rayleigh 
calculated the most unstable wave to have wavelength approximately nine times the jet radius. 

ven though breakup into drops is a nonlinear phenomenon linear theory does surprisingly well 
m predicting overall features such as drop sizes and breakup times. For a quantitative description 

of breakup linear theory must be abandoned and we shall be exclusively concerned with nonlinear 
solutions here. 


This classical problem has received considerable attention over the years. Experimental work 
has been carried out by Donnelly and Glaberson [10], Goedde and Yuen [15] and more recently 
by Chaudhary and Maxworthy [7], [8], The related problem of the dripping faucet was studied by 
Hauser et al. [16] and in much more detail by the recent work of Peregrine, Shoker and Symon 
IJIJ and [40J, [3] where sequences of photographs detailing the breakup can be found. The latter 
authors report on a fascinating cascading phenomenon during pinching. 

Theoretical studies including nonlinear effects are more recent and roughly fall into three cate- 
gories: Direct numerical simulations, weakly nonlinear theories and slender jet theories Since we 
will be concerned with both inviscid and viscous flows the literature is reviewed separately. 

For inviscid jets, theories stemming from the classical weakly nonlinear approach of Stuart [42] 
and Watson [46], have been applied by Yuen [48], Chaudhary and Redekopp [9] and Lafrance [211. 
These approaches follow the evolution of two low amplitude initial perturbations, a fundamental 
an a harmonic, and an amplitude equation with cubic nonlinearity is derived in the classical way. 
An objection to this methodology for the jet problem is that linear stability always gives a band 
of unstable waves which is unchanged as physical parameters are varied, in this case the capillary 
number. The evolution equation, then, is only valid for small enough times and in particular 
cannot describe pinching since that requires amplitude levels beyond the theory’s range of validity. 
Chaudhary and Redekopp used their amplitude equations beyond their range of validity to produce 
pictures that qualitatively agree with experimental observations of breakup; there is no reason 
however, to expect any quantitative agreement and in particular accurate solutions near breakup ’ 
An improvement from weakly nonlinear theories is to allow the interfacial amplitude to be 
arbitrary but its wavelength to be long compared to the undisturbed jet radius This is in the 
spirit of shallow water wave theories (see Whitham [47] for example), and through a systematic 
asymptotic expansion one-dimensional models of the Euler equations arise; the small parameter 
used is the ratio of the undisturbed jet radius to a typical perturbation wavelength. This was 
done by Ting and Keller [44], referred to as TK, who derive and study the leading order system. 
TK construct similarity solutions of their equations and use them to study the flow just beyond 
pinching by improving the original arguments of Taylor [43] and Keller [18]. A further improvement 
has been recently added by Keller, King and Ting [19] who solve the Euler equations for the flow 
in the attached fluid blob and match it to the slender filament which is a solution of the slender jet 
equations. It has been found out recently (see Papageorgiou [29], Papageorgiou and Orellana [30]) 
that the initial value problem based on the TK system terminates in an infinite slope singularity 
before the jet radius vanishes to form a pinch (see later also). It seems unlikely, therefore, that 
the leading order equations can describe pinching and keeping within the context of long-waves 
additional terms in the expansion need to be kept. Such models have been known for some years 
and are collectively called Cosserat models. A review can be found in Bogy [2], Meseguer [25] and 
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a more recent work with additional references is that of Garcia and Castellanos [14]. It appears 
that one of the original computations based on a one-dimensional model is that of Lee [22]; it was 
found that the Lee model supports pinching solutions but as our numerical studies indicate the 
jet develops a cusp near pinching - this in turn implies the violation of the slender jet assumption 
at the time of pinching and the full Euler equations need to be incorporated in an inner region. 
A weakly nonlinear version of the Cosserat model has been used by Schulkes [35] as a basis for a 
computational study. Schulkes [36] makes an evaluation of the weakly nonlinear model with the 
Lee model by solving the one-dimensional systems numerically; the local behavior very close to 
pinching was not studied, however. 

Even though the one-dimensional models are useful and desirable, it appears that for inviscid 
jets at least, it may be necessary to relax the slenderness assumption and work with the full 
Euler equations. This is the axisymmetric extension of the work of Keller and Miksis [20] who 
compute the evolution of a break in a fluid sheet (the slender theory of this can be found in 
TK). A similarity solution of the Euler equations with all the terms balancing as a pinch forms, 
is possible then but has not been computed yet - see Papageorgiou [29]. Such solutions would 
be useful in comparisons with simulations based on the Euler equations. Such computations are 
usually based on boundary integral methods and have been performed by Mansour and Lundgren 
[24] for breakup of a liquid column and for related problems such as pendant drop bifurcations 
and contraction of liquid filaments by Schulkes [37], [38], [39]. A detailed quantitative comparison 
between computations and local similarity pinching solutions has not been carried out yet. 

The difficulties encountered in inviscid one- dimensional models are removed when viscosity is 
added in the sense that the slenderness assumption is valid all the way to pinching governed by 
the leading order equations. Both Stokes and Navier-Stokes jets are considered. In the absence 
of inertia a one- dimensional model of the Stokes equations was first derived by Renardy [32] who 
also includes various viscoelastic models. Using the slender jet equations, Renardy proved that 
a Newtonian jet will break after a finite time given appropriate initial conditions. The solutions 
near the singular time are self-similar and have been given by Papageorgiou [28] who obtains an 
implicitly defined closed form solution with universal scaling exponents; in addition the scaling 
functions are universal to within a stretching factor that comes from initial conditions. All the 
theoretical results were confirmed by numerical solutions with different initial conditions. The 
jet shape is locally symmetric near the pinch point and the minimum jet radius goes to zero like 
(f s _ t) 0 with /3 « 0.175 and t s the singular time which depends on initial conditions. Other 
forms of breakup similarity solutions are given by Renardy [33] and correspond to different initial 
conditions. In the study of Papageorgiou [28] the initial conditions are smooth. 

Direct numerical simulations based on the Stokes equations and the boundary integral technique 
are usually carried out by inclusion of a surrounding fluid. In fact the single jet problem is usually 
computed by taking the appropriate viscosity ratio limit. Such computations have been developed 
by, for example, Stone and Leal [41] and Tjahjadi, Stone and Ottino [45]. An evaluation of the self- 
similar solutions of Renardy and Papageorgiou using the direct simulations remains to be carried 
out. 

Inclusion of inertia into the Stokes theories enables construction of pinching solutions based on 
the Navier-Stokes equations. The analogous set of evolution equations was first obtained by Eggers 
and Dupont [13] who used local Taylor expansions assuming slenderness. An equivalent derivation 
has been given by Papageorgiou [28] using a slight modification of the systematic asymptotic 
expansion procedure applied to Stokes jets, by increasing the Reynolds number away from zero 
and scaling it with the expansion parameter. It should be noted that the model of Eggers and 
Dupont includes the complete curvature terms in the normal stress balance even though not all 
terms there are of the same asymptotic order. Eggers [11], [12], has computed the self-similar 
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solutions near breakup according to the one-dimensional model and provides strong numerical 
evidence that the scaling exponents and scaling functions are universal. In addition the breakup is 
non-symmetric in general if the initial conditions are also non-symmetric - this is in contrast to the 
Stokes models which produce locally symmetric pinching jet shapes even if the initial conditions are 
non-symmetric. It appears that the one-dimensional models do very well in describing the pinching 
dynamics and could be used in hybrid numerical methods to continue computations beyond a 
change in topology. It would be interesting to compare the theoretical predictions with direct 
numerical simulations of the Navier-Stokes equations. There have been some recent advances in 
this direction and we cite the work of Richards, LenhofF and Beris [34], and a review can be found 
in Beris et al. [1], 

In the theories described above the analytical approach is divided into two steps: First, a 
simplified set of one-dimensional PDE’s is derived by assuming slenderness; second, these PDE’s are 
analyzed for pinching self-similar solutions and determination of scaling exponents and functions. 
In solving initial value problems, then, the jet is assumed to be slender throughout the evolution. It 
is possible, however, that the evolution prior to pinching does not obey the slenderness criterion but 
only does so near pinching and locally near the pinch point. One can look for local solutions of the 
full equations of motion with a slenderness ratio which is time dependent and in fact approaches 
zero as the singular time is reached. The two steps are combined into one and the final result 
is a set of ODE’s for the self-similar scaling functions; these are of course the same as the two 
step approach. This method avoids any problems with transient behavior of the evolution systems 
and gives a clearer structure of the singular solutions of the full equations. Such constructions 
starting from the Euler, Stokes and Navier-Stokes equations can be found in Papageorgiou [291 and 
is described later. 

The structure of the article is as follows. Section 2 is devoted to the study of breakup according 
to the Euler equations. One- dimensional models are considered first and evidence is given suggesting 
the need of non-slender inner solutions; the equations governing such inner solutions are also derived 
and partly analyzed. Section 3 looks at Stokes jets and includes the analysis of pinching solutions 
both from the one-dimensional model as well as directly from the equations. Numerical solutions 
are presented that confirm the theoretical predictions. A modification of the expansion in Section 
3 leads to jets with inertia and a model identical to that of Eggers [11]. We show how the universal 
scaling functions can be used to continue the dynamics beyond pinching. This is included in Section 
4. Section 5 includes conclusions and thoughts on future work. 


2 Inviscid jets governed by the Euler equations 


Using a cylindrical coordinate system and looking at irrotational axisymmetric flows, the equations 
of motion can be written in terms of a velocity potential (j>(t,r,z) where t is time, r is the radial 
and z the axial coordinate. The corresponding velocity field is u = (u, v, w) = \7<f>. The governing 
equations and boundary conditions are 


On r = S(t, z) 


<f>TT H <j>T + 4>ZZ = 0. 

T 


(i) 


<t>r 

4>t + + <t>l) 


St + 4>zS z , 



s 2 z s„ 
1 + 52 
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(1 + S Z 2 ) 1/2 . 


( 2 ) 

( 3 ) 



In addition we require regularity of velocities at r = 0. The equations and boundary conditions 
have been non- dimension alized using the undisturbed jet radius R as a length scale, the quantity 
(pU 3 /ct) 1 / 2 as a time scale (this is the time scale of capillary instability), while the potential is scaled 
by ( Rc/p ) 1 / 2 ; p and a denote fluid density and surface tension coefficient between the liquid and 
air. If initial and boundary conditions are prescribed (for instance periodic boundary conditions) 
the system (l)-(3) is closed and must be addressed numerically in general. In what follows we 
consider possible analytical developments. 

2.1 Slender jet theory: One-dimensional models 

Consider a small positive parameter c which physically measures the ratio of a typical interfacial 
amplitude to a typical wavelength. The following change of variables gives the long wave equations 
to be studied, 

d r — ► d r , d z -> ed z , d t — *■ ed t , 4> -» -<t>, 

which casts (l)-(3) into 

<j>TT H 4>T + ^4 > ZZ-> (4) 

r 

on r = 5(/,z) 

4>r = e 2 (S t + 4>zS z ), (5) 

* + s($ + *) - (6) 

An asymptotic expansion for <f> proceeds in powers of e 2 

<f> — <t> o + £ 2 4>\ + • • • i 

and the leading order solutions follow from (4): 

4>0 = 4*oi4i 4*1 = 4*0zz- ( 7 ) 

Using these solutions into (5)-(6) yields the following system correct to 0(e 2 ): 

St + —S4>o zz + S z 4*0z = 0 , ( 8 ) 

4*0t + 2 ^ 0 z + ~ ^ S 2 (4>q zz — 24>Oz4*Ozzz ~ %4*0zzt) = ~~g + f2 Szz + € 2 ^. • (9) 

Next, defining W = <p Qz and differentiating (9) by z we get the following system governing the 
evolution of the jet shape and the axial velocity: 

(5% + {S 2 W) Z = 0, (10) 

W t + WW z + U 2 (s 2 (W* -2WW zz -2W zt j) z = |§ + e 2 S zzz + ^e 2 (^J . (11) 

For a more detailed derivation as well as the inclusion of a surrounding fluid see Papageorgiou and 
Orellana (1995). 
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When € = 0 equations (10)-(11) are equivalent to the system studied by TK. It has been found 
by us, however, that this system terminates in “shock” finite time singularities and consequently the 
slender jet assumption is violated before the radius vanishes. The evidence is numerical by solution 
of the initial value problem on periodic domains using spectral methods, and also analytical by 
treating the system as a 2 x 2 system of conservation laws which becomes hyperbolic in the complex 
plane when the dependent and independent variables, but not time, are complexified. Riemann 
invariants are used to show that there are envelopes of characteristics which reach the real axis 
and thus a shock forms after a finite time. For small amplitude initial conditions the numerically 
computed singular times are in excellent agreement with the theoretical predictions. For full details 
of this analysis and the accompanying computations see Papageorgiou and Orellana (1995). Much 
of this theory stems from the pioneering work of Moore [26], [27] on singularity formation in vortex 
sheets and the subsequent work of Caflisch and Orellana [5], [6], Also of interest is the related 
recent work of Caflisch et al. [4]. 

Since our interest is in drop formation we contemplate next regularizations of the TK equations 
as in (10) and (11). In particular we will present here the analysis of the Lee model which serves as 
a good illustrative starting point in the evaluation of one- dimensional inviscid models. We note that 
this model has been used by several investigators including Lee [22], Meseguer [25] and Schulkes 
[35], [36] - see Introduction also. The Lee model derives by dropping the c 2 terms from the left 
hand side of (11) but keeping the full curvature terms on the right. The model considered here 
does not keep the full curvature terms but ignores the t 2 terms on the left. Preliminary numerical 
results show that the retainment of the full curvature term has little effect on the final structure 
of the singularity. The model to be studied next, which is also known as a slice model, is 

(S 2 )t + {S 2 W) Z = 0, (12) 

W.+WW, = + . (13) 

We note that without loss of generality we can scale c out of the problem by the transformations 

S -* -5, W-^^eW, z —* t. 

€ y € 

Mathematically, these transformations do not alter the form of the singularities which we are 
analyzing. 

2.2 Finite time singularities 

Assuming that after a time t s the jet radius vanishes at some position z $ , we look for solutions near 
this time and position, of the form 

S(z,*) = r“h(£), W{z,t) = T‘ 3 -'g{t), (14) 

where r = t s — t > 0. We note that the scaling for W arises from a balance of terms, as r — 0, 
in the mass balance equation (12). The value of a depends on which terms are in balance in the 
Bernoulli equation (13) near the singular time. A reasonable assumption would be to assume that 
all terms in (12) are in balance as r -»• 0 which in turn fixes that a = (i = 2/3. Similarity solutions 
can then be constructed by study of the corresponding ODE’s for the scaling functions. We have 
found, however, by numerical solution of the evolution equations (12) and (13) that the most likely 
scenario is one that has 


(5 > a. 


(15) 
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Applying this inequality and an orders-of-magnitude argument to (13) it is found that the term 
S 2 /S 2 on the right is negligible as r — ► 0, and the balance gives 

a = 4/3 -2. (16) 

Substitution of the ansatz (14) into the evolution equations along with the inequality (15) yields 
the following equations for the scaling functions 

—4(2/? - 1 )h 2 + 2 fchti + h 2 g’ + 2 ghti = 0, (17) 

(1 - 0)g 4- (/?£ + g)g' - h!" + - . (18) 

The inequality (15) implies that as the jet radius vanishes the interface develops a cusp there 
and so the slender jet ansatz is violated. The situation is better than in the TK model, however, 
since that model does not support the vanishing of the jet radius before the long wave assumption 
is violated. It is routine to compute representative scaling functions - a four parameter family of 
similarity solutions is obtained by specification of p(0),h(0),h'(0),h // (0) and integration to ±oo. 
A subset of this family has symmetry in h and anti-symmetry in <7, which is a property of the 
equations. For this case, the asymptotic behavior at infinity is easily found to be 

MO ~ A^- 2 /3, g{ o-r 1 ^, 

with A and B constants. Next order corrections can be calculated and decaying oscillations at 
infinity are found by application of a WKB method - see [30] for details. /.From the numerical 
solution of the evolution equations (see next subsection) it is found that the most likely case is 
/? = 0.6, a = 0.4. These values imply h ~ £ 2 ^ 3 an d this power law was confirmed with solutions 
of (17)-(18) which give the exponent value to be 0.668. (The conditions for these solutions are 
g(0) = 0, h( 0) = 1, h'( 0) = 0, h"{ 0) = 0.5.) 

2.3 Numerical solution of the evolution equations 

Equations (12)-(13) were solved numerically using accurate pseudospectral methods (see also [30]); 
as explained earlier, we take e = 1. Periodic boundary conditions are imposed and in addition 
(this is not built into the numerical scheme) the symmetries of the equations are utilized by giving 
initial conditions S(z, 0) = 5o(z) and W(z, 0) = Wq(z) which are even and odd respectively about 
z — tt where the spatial domain is [0,27r]. This implies that S(z,t ) and W(z,t ) remain even and 
odd respectively for subsequent times. The initial conditions used in the results that follow are 

S 0 (z ) = a + ai cos (z), W 0 (z) = -fesin(z), (19) 

where a > |aj|. The integrals Jq* W(z,t)dz and S 2 (z,t)dz are conserved quantities and were 
used to monitor the accuracy of the computations. 

The evolution of a typical run having b = 0.1, a = 0.5 and ai = 0 is shown in Figure 1. It 
is seen that drop formation appears after a finite time with the mini m um jet radius going to zero 
at two symmetrically placed positions with the corresponding axial velocity blowing up there, as 
expected. (Fluid is ejected out of the drop which is forming.) The computations were stopped 
when (i) the minimum jet radius dropped below 0.01, (ii) the value of max|W| exceeded 20, or, (iii) 
the value of max|Wz| exceeded 100. Improvements can be made by inclusion of more modes but 
it is found that these criteria are sufficient to infer singular times and scaling laws, for example, 
described next. 
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2.4 Singular times and scaling laws 

The time of the singularity is an unknown quantity and depends on the initial conditions It is 

7 I 1 " 6 " ae,1 ' 0d ° f **■“»• “ estimate of this itt orde7o use I, 

. • 4 1 ° f SCabn§ ex P° nents - The following method appears to be the most suitable one 

and in particular it is insensitive to the type of ansatz we wish to evaluate. From the analysis 
presented earlier we expect that as a pinch forms the value of W z blows up as r" 1 independently 
of the values of the exponents a and 0. Letting max|W z | = W m (t) we have 

WjJ) ~ ^ ~ ^ 35 t — ► t s — . (20) 

Plot « of ^ produced linear behavior by the end of the computation and a least squares fit of 
this hnear part was used to extrapolate and estimate This procedure is depicted in Figure 2 
which corresponds to initial conditions with a = 0.5, a a = 0 and 6 = 0.05. Singular times for 
different initial conditions with 6 varying were also computed. 

<14vTr i^M a V« meS c n0 ^ ‘ Um MXt *° tie es ‘ imatio11 "f the exponents o and 0 (see 
(14)). It is useful to define S„ m (t) = min(S(z,t)) and W„.„(i) = max(W(z,t)), noting that these 

“ S r r ° n ' , P '° ,S ° f ^ ‘"(H'J) against in(r) fhould yield 

straight hnes near the singular time (i.e., as r - 0) with slopes equal to a and 0 - 1 respectively. 

igures 3a- b depict such estimates for the representative case b = 0.25. Clearly a < 0 and so a 

cusp singidarity is forming. As described earlier a balance of terms in the Bernoulli equation as 

ofTnitil! v T $ P ' ( 16 > betwee * « ^d 0. It is important to test this balance for a range 

a are 05^“^ nil' f f ° r ^ eXample ° f FigUre 3 ’ the estimated values of 40 - 2 and 

summarv f” °' 44res P ectlvel y and are fairly close given the sensitivity of the computations. A 
summary of results for various initial conditions is given in Table 1 below which contains estimates 

bda ’ ’(? 6 )“ aS ^ ~ WhKh Sh ° uld be com Pared with the value of a as indicated by the 


b 

.005 

.01 

.025 

.05 

.075 

.1 

.125 

.15 

.175 

.2 

.225 

.25 

.275 

.3 

u 

5.89 

5.10 

4.05 

3.25 

2.73 

2.29 

1.94 

1.66 

1.44 

1.25 

1.08 

.938 

.810 

.695 

a 

/j 

.47 

.48 

r.47 

.46 

.46 

.46 

.49 

.47 

.47 

.47 

.45 

.44 

.43 

.43 

P 

A O O 

.61 

.53 

.61 

.59 

.59 

.61 

.64 

.61 

.6 

.62 

.6 

.63 

.56 

.57 

4p — 2 

nni 

.44 
1 1 • 

.12 

.44 

.36 

.36 

.44 

.56 

.44 

.4 

_ .48 

.4 

.52 

.24 

.28 


ship * tuiiiMuuu ana mey are consistent with the relation- 

Other features of the dynamics can be quantified using the model equations. For instance it 
has been shown that breakup according to the most amplified linear wave (this is an empirical 
engineering estimate) underestimates the singular time. More interestingly, the model predicts 
drop formation and can be used to compute the drop volume for different initial conditions. It has 

been found (for the symmetric initial conditions used here at least) that the drop volume roughly 
vanes according to the estimate ® * 

Vdrop — 0.51 exp( — 8.66). 

This is estimated from our computations and interpolation of the data. 

2.5 Non-slender breakup 

The cusp singularities at breakup violate the slender jet ansatz and an inner solution must be sought 
which does not assume scale disparity. Accordingly, the following transformations provide exact 
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similarity solutions of the Euler equations in the sense that all terms in the governing equations 
aTe in balance then (note that such an ansatz follows from a dimensional analysis also, as applied 
by Keller and Miksis [20] for the breakup of a fluid sheet). 

(r,z) = r 2/3 {y, 0, S(z,t) = r 2/3 /(0> 4>(r,z,t) = T^ 3 X (y,0- ( 21 ) 

Substitution of (21) into the governing equations (l)-(3) yields the following problem to be solved: 

Xyy + ~Xy + X« = 0. ( 22 ) 


On y — /(*), 


l(2ZXi + 2fXv-x) + kxl + Xi) = ~ 


Xy = §(«'-/) + /'%. 

1 f „, /*/" 
/ 1 + 1+/' 2 


(i + /'T 1/2 - 


0\— /vy / 2' 

Note that equations (22)-(24) are the axially symmetric analogue of those used in [20]. 


(23) 

(24) 


3 Viscous jets 

When viscosity is important the pertinent scales for capillary instability are <r//x for velocities, a/R 
for pressure and pR/cr for time, where p is the kinematic viscosity of the fluid. The non-dimensional 
form of the radially symmetric Navier- Stokes equations together with boundary conditions are: 


u 

R e (ut + uu T + uw z ) = -p T + Ati - 

R e (w t + uw r + ww z ) = — p z + Aw, 

u 

u T + - + w z = 0, 

r 

where A denotes the Laplacian in cylindrical coordinates. On the interface r = S(z,t): 


(■ U Z + uv)(l - Si) + 2 u r Sz - 2 WzSz = 0 , 


p(l + Si) - 2 u T - 2w z S z + 2(u z + w T )S z = - 




L\A + si s 

u = St + wS z . 


(25) 

(26) 
(27) 


(28) 

(29) 

(30) 


In addition to the interfacial conditions (28)-(30), we require boundedness of solution at r = 0. 
The non-dimensional group R e = ( cpR/p 2 ) is the reciprocal of the capillary number, C a . 

In general a numerical solution of this system is required and possible approaches, when a 
surrounding fluid is present, are described elsewhere in [34] and also the article of Beris et al. in 
this volume. In what follows we use asymptotic methods to tackle the problem of jet pinching, 
viewing such an event as a finite time singularity of the PDEs. Such analyses are of value in 
evaluating direct numerical simulation and in the continuation of the dynamics beyond pinching 
by rational rather than ad hoc approaches. 
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3.1 Slender jet theory 

The ideas foUow those of the inviscid analysis where essentially we transform the governing equa- 
tions according to d/dz — ► ed/dz. In order to bring in unsteady and nonlinear terms into the 
leading order evolution the following ordering of R e is required, 

R e = € 2 k, k = 0(1). (31) 

When k = 0 we recover the slender jet equations for Stokes flow. An asymptotic solution is sought 
in powers of e 2 as e — *• 0 in the form 


u(t, r, z) = u 0 + e 2 u a + . . . , 
w(t, r, z) = -w 0 + (Wi + . . . , 

(32) 

(33) 

p{t,r,z) = po + c 2 pi + . . ., 

(34) 

S(t, z) = 50 + i 2 S\ + . . . . 

(35) 


Substituting (31)-(35) into the governing equations (25)-(27) and collecting terms gives a series of 
simpler problems to solve. To leading order the axial momentum equation (26) and the continuity 
equation (27) provide the leading order flow field in terms of a single unknown function of z and t: 

w 0 = w 0 (t,z ), u 0 = -^rw 0z . (36) 

At the next order equation (26) gives 


K(w 0i + U> 0 *%j:) = ~POz + W lrr + ^U) lr + W 0zz , (37) 

while at leading order the radial momentum equation (25) brings in the radial variations in pressure: 

POt — u Orr + ~ u 0 r (38) 

Substitution of the solution for uo (36) into the equation for the leading order radial pressure 
gradient (38) shows that the latter is zero and so is independent of r. We write, then, 

Po =po(t,z). 

Since the pressure does not vary radially, to leading order at least, it must be equal to its value 
just below the jet surface. The normal stress balance (29) allows its evaluation (using the leading 
order flow field already calculated): ' 


Po = j- o ~w 0z . (39) 

We have two unknowns, So and wo, and have yet to use the tangential stress balance (28) and the 
kinematic condition (29). A system of leading order evolution equations is found as follows: We 
first solve for the axial perturbation velocity uq from equation (37); this is easily done since the 
forcing terms involve w 0 and p 0 which are both independent of r. With uq at our disposal we can 
evaluate the two leading orders in the tangential stress balance equation (28). The leading order 
contribution is trivial ( Wo r = 0) while the next order contribution gives a PDE involving wq and 
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So- Secondly, the kinematic condition is known to leading order since the flow field is also known 
to that order. The following coupled system of PDEs needs to be addressed, then: 


k(wq t + WqWoz) = 


3(SqWoz)z 

sl 



Sot + + w oSqz — 0 . 


(40) 

(41) 


The system (40)-(41) must be addressed numerically in general. We note that k can be scaled out 
of the problem but is retained in order to clarify the Stokes limit considered next. A derivation and 
a numerical study of (40)-(41) was first carried out by Eggers and Dupont [13] using a different 
but related expansion scheme and additional work can be found in Eggers [11], [12], Papageorgiou 
[28], [29] and Renardy [32], [33]. The latter author considers inertialess flows. In what follows we 
consider the two distinct limi ts of Stokes flows (no inertia k = 0) and flows with inertia. 


4 The model equations for Stokes flows 


With k = 0 equation (40) can be integrated once to give 


Wo * = l 


rm 

Uo 2 



(42) 


where the function A (t) represents the force in the jet (this force is independent of z since the 
flow is inertialess - see Renardy [32] and Papageorgiou [28]) and is to be found. If in addition we 
impose spatial periodicity, as will be done in the numerical experiments presented later, a further 
integration of (42) gives the dependence of the quasi-uniform force in terms of the jet shape, 


_ S?(USo)dz 

U JT(1 !Sl)dz' 


(43) 


where, without loss of generality, the period has been normalized to 27r. This expression is useful 
in the numerical work given later. 


4.1 Finite time singularities and self-similar solutions 

Assume that the jet pinches after a finite time t a . Then the minimum jet radius vanishes at t = t a at 
some axial position z = z*. The objective is to obtain a description of the flow in the vicinity of this 
space-time singularity. Denoting t = t s — t, we construct solutions in the limit 0 < r < 1 near the 
position z = z a , i.e. as \z-z a \ — ♦ 0. Working with (40)-(41) and employing an order-of-magnitudes 
estimate (see Papageorgiou [28]) we may look for similarity solutions of the form 

S 0 (t, z) = r/(£), w 0 (t, z ) = r^ -1 5 (f), f = (z - z a y e (44) 

where 0 < 0 < 1 and the scaling functions /(£) and g(£) are to be found. Substitution of (44) into 
(40) and (41) with k = 0, yields the following set of nonlinear ODEs to solve on f G (— oo, oo): 

4(S/V + /) = 0, (45) 

VS 

(ff + /W' + (^-l)/ = 0. (46) 
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As before, equation (45) may be integrated once to obtain 


, 1 k 

9 3f + p* (47) 

where k is a constant of integration which is expressible in terms of / once we note that g(£) — > 0 
as |f | — ► oo; integration of (47) gives 


, _ i FLWfW 
3/!L(i/m‘ 


(48) 


This expression for k can be used to obtain the form of X(t) as t t s . Substitution of the similarity 
scalings (44) into (42) and integration of the resulting expression with respect to £ using the decay 
of g at infinity, leads to 


A = 3 k(t s — t) as t — ► t s — . (49) 

This assertion follows from the analysis but needs to be checked from a solution of the initial value 
problem for (40)-(41) for k = 0. 


4.2 Solution of the similarity equations 

It is easy to establish the behavior of / and g for large |£|. This can be done by directly balancing 
terms in the ODEs or equivalently noting that the solutions at infinity be independent of t when 
written in terms of outer variables. Both approaches lead to 

/(0~I£| 1//J , 9(0 ~ \t\- {1 - P)/l3 as |£|-oo. (50) 

In terms of outer variables this asymptotic behavior provides matching conditions So r*j 1 2, — Z s |i/j0 
and wo ~ \z — z s as | z — z s \ — 0, the multiplicative constants coming from solution of 
the system (45) and (46). Since g(£) tends to zero at infinity, there is a point, say, where 
g(£o) + flf,o = 0. Applying the transformations 


/ — /, <3(7?) = g + (3£ o, r? = ^ - £ 0 , 

shifts this point to the origin, r? = 0, where (3(0) = 0. The transformed equation (46) has a singular 
point at the origin which is removable if <3'(0) = 2; the solutions are smooth, then. A local analysis 
of the equations for small 7? gives 

/(*?) = /o + t ? 2 /2 + t? 4 /4 + • • • , Gfa) = 2 t? + r) 3 g 3 + . . . , (51) 


where we find 


fo = 


1 

12(1 + py 


3 + 2/3 

72(1 +/3) 2 ’ 


(52) 


with the remaining coefficients /, and <?,-, i > 2, expressible in terms of the single parameter f 2 . 

An implicitly defined closed form solution is available after elimination of (3 to obtain a second 
order ODE for /. This solution is (for details see Papageorgiou [28]) 


1 

[12(1 + /3)F 


l 


f = i 

J 12(1 + (3) 

6 ( cosh 2 6 + 3 + 2/3)^ +1 /2 


cos h 2 (9). 


cosh^ 


d9 — ±Ar], 


aMm L(-Tf + 7)* ,+Mk 


(53) 
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where 


A = 


24 f 2 + 13 ^ +1/2 T 
(l + (3)\6(l + /3)J V/2 ’ 


and ± corresponds to 77 positive and negative respectively. Evaluation of the solution G{rf) above, at 
±00 and use of the boundary conditions G(-oo) = G(oo) = /?fo gives an expression for k identical 
to (48) noting the symmetry of /. The ratio is independent of /2, so given (3 the solutions / and 
g are obtainable numerically. An eigenrelation must be satisfied, however, in order to obtain the 
correct (3 and the corresponding scaling functions. This relation arises from the two expressions of 
k given by (48) and (52) and reads 

3+2/? _ 1 

72(1 + /!)* 3/ 0 °°(l//W 10 ' 

The computations give (3 = 0.175 correct to three decimals and this value appears to be unique. The 
effect of /2 is in rescaling the axial coordinate as can be seen from the sample numerical results 
shown in Figure 4. The local theory provides similarity solutions which are unique to within a 
linear scaling transformation. The theory is a formal one and numerical solutions of the initial 
value problem are desirable to confirm the local theory. This is done next. 


4.3 Numerical solutions of the initial value problem 

The numerical methods employed to solve (40) and (41) for k = 0 are described in Papageorgiou 
[28]. A pseudo-spectral method is used with axial periodicity; the boundary conditions are insensi- 
tive to the terminal state solutions since the latter are indeed of the local form given in the previous 
subsection. To fix matters we consider periodic solutions on [0,27r]; the equations preserve sym- 
metry and anti-symmetry respectively for So and wo about z — tt, and in the first set of numerical 
results we report on such parity solutions. The following initial conditions are used: 

5 o (0, z) = a + b cos(z), (55) 

where a, b > 0 and a > b. This means that there is a minim um in the jet shape at z = 7r and due 
to symmetry the jet will first pinch there. In addition to So and wo we also monitor the evolution 
of A(£) and the minimum jet radius denoted by 5 m , n (£). A typical run having a = 0.5 and 6 = 0.1 
is given in Figure 5 which shows evolution of the jet shape and the corresponding axial velocity at 
different times up to t = 6.5. A cross-section of the jet is shown at the later time t = 6.6 in Figure 
6; the Figure also shows the evolution of the maximum value of So* and confirms the validity of 
the slender jet assumption. Note that the singular time (see below for its estimation) is t $ = 6.64 
and the minimum jet radius at t = 6.6 is 0.0028. 

According to the predictions of the local theory, as the singular time is approached the following 
forms are attained by S mtn (f) and A (t): 

1 q _l 0/3 

s ^~wr «('■-')' x u~whw^- tl (56) 

where (3 = 0.175. Considering either of these two expressions, then, we can estimate the singular 
time t s since the behavior is linear near the singular time. This procedure is illustrated in Figure 7 
using data from the previous run. The Figure shows the evolution of S m i n (t) and an enlargement 
of the later stages of the evolution together with the least squares linear extrapolation used in 
estimating t s . This is how the value t s = 6.64 was found. It is worth noting that the variation 
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°f S m i n (t) is linear long before the singular time, i.e., the asymptotic behavior predicted by (56) 
for small t s — t is valid for most of the computational time - in fact from approximately t = 2 to 
t = t s . Similar findings hold for a parallel analysis of the data from the X(t) evolution (see [28]). A 
further confirmation of the wide validity of the asymptotic solutions and the universal exponents, 
is given in Figure 8 which depicts the computed and theoretical evolutions of S min (t) and X(t); this 
is possible once t s has been estimated and for the theoretical curves we use the results (56). Again 
agreement is seen to be excellent with the similarity solutions valid significantly before breakup. 

The corresponding scaling functions are computed as follows: Knowledge of the singular time 
t s , and hence r = t s - t, allows the computation of the similarity independent variable f as the 
computation proceeds - this is given by ( z-ir)/r * for the run described earlier. The corresponding 
computed solutions at that particular time, So(t, z) and wo(t, z), are next used together with the 
similarity ansatz (44) to construct the functions /(£) and g(£). The dependence of these functions 
on r disappears as r — ► 0 and this expectation is fully supported in the sample numerical results 
of Figure 9 which shows the corresponding /(£) and p(£) for different times approaching t s . 

Having established the universal predictions of the similarity theory for solutions containing 
special symmetry, we next consider breakup from more general initial conditions. In particular we 
report on a run having 


<5b(0, z) = 0.5 + 0.1(sin x + cos x). 

The evolution of So and Wq is given in Figure 10. The jet is breaking at a position z = z s (t) to 
the right of z = x. An estimate of the singular time can be made as in the symmetric case (the 
similarity theory is a local one) and the evolution of S min (t ) and A(f) together with the theoretical 
predictions are shown in Figure 11. Again agreement is excellent. We also note that the breakup 
time is smaller. The theory presented in the previous subsection predicts scaling functions /(£) 
and $(£) which are even and odd respectively about f = 0. This prediction has been confirmed 
for this nonsymmetric run as follows: Take the last computed profiles So and wq. From the 5b 
profile locate the point where the minimum is attained and call this z 0 . Then S 0 is expected to be 
symmetric about z = zo and wq anti-symmetric. This is presented graphically in Figure 12 for So 
and wo respectively. The Figure for So was obtained by reflecting the solution for z < zq about 
z = zo and plotting it along with the solution for z > z 0 , while that for w 0 is obtained by reflecting 
the solution about z = z 0 , taking its negative, and plotting it together with the solution for z > z 0 . 
Solid lines represent z > z 0 and the open circles correspond to reflected data. The two data sets 
are indistinguishable providing strong evidence for the local symmetry properties predicted by the 
similarity theory. 

4.4 Universal scaling functions 

In the previous subsections we established a way of (i) estimating the singular time t s , (ii) using 
this estimate and numerical solutions near breakup to construct the scaling functions given by (44). 
In this section we consider the scaling functions for different initial conditions. Since breakup is 
locally symmetric we consider a one-parameter family of symmetric initial conditions of the form: 

5o(0, z) = 0.5 + b cos z, 

Runs with b = 0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3 are reported. Figure 13 shows the evolution 
°f for the different values of b labeled on the Figure. In addition the asymptotic behavior 

of S m in near the singular time is superimposed for each run; the slopes of the curves are predicted 
by the theory to be —1/12(1 + 0) irrespective of initial conditions, and this is supported by the 
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results. Another interpretation of the equal slopes is that the /(£) scaling functions satisfy /( 0) = 
1/12(1 + /?) for the different values of b. 

A more detailed study of the scaling functions is considered next. Using the methods of Section 
4.3, scaling functions were constructed for different values of b. Representative results for /(£) are 
shown in Figure 14a for b (denoted by ep in the Figure) equal to 0.005,0.1,0.2,0.3. The scaling 
functions coincide at £ = 0; according to the theory of Section 4.2 and in particular the solutions 
(53), the only difference between scaling functions is in a re-scaling of the axial coordinate brought 
about by the presence of the constant / 2 which depends on the initial conditions. This is illustrated 
numerically from the data of Figure 14a as follows: Consider the run with b = 0.1 as a reference run 
and call the emerging scaling function fo(r)). Then, according to the theory, there exists a positive 
number c(6) such that the transformation r) —> cr) maps the scaling function /(£; 6) into / 0 (£). The 
value of c can be obtained by consideration of a fixed interfacial height, h say, and is equal to the 
ratio 770 / 77 ,' where fo(vo) = = h with i corresponding to the different initial conditions. The c 

found from mapping a single point onto the curve / 0 (»7) will map the whole curve onto /o( 77 ), if the 
theory and the numerics coincide, and the procedure outlined has the dual purpose of (i) making a 
strong global validation of the similarity solutions, (ii) constructing “universal” scaling functions, 
noting that these are so to within a linear axial transformation. The results of this procedure are 
given in Figure 14b from which it is seen that all curves do collapse into the single fo(ri) curve, in 
complete accord with the theory. 

4.5 Pinching solutions of the Stokes equations 

The Stokes equations for axially symmetric liquid threads are given by (25)-(27) with R e = 0 and 
the same boundary conditions (28)-(30). Here we show how the two step procedure outlined above, 
that is the derivation of longwave equations followed by construction of pinching solutions, can 
be combined into one by essentially allowing the small parameter e introduced into the longwave 
models, to become time dependent. Using the notation introduced earlier, we generalize the breakup 
by looking at solutions of the Stokes equations near the singular time, i.e., as r — ► 0, of the form 

r = r a y , z = A, S = T a f(0, 

w = rW(f,y,0, *sr 7+a "^l7(«,»,f), P = r _ “P(f,t/,£), (57) 

where without loss of generality we take the breaking point to be at z = 0. Note that the form of 
u follows from the continuity equation once w is assumed. We make the assumption that a > (3 
which implies a slender geometry which is enhanced as the breaking time is approached; slenderness 
is not assumed throughout the evolution, then, but only near the pinch time. The quantity r 2a-2/3 
is a small parameter, then, which appears naturally in the Stokes equations (in the biharmonic 
operator for the streamfunction, for instance). The following Taylor expansions are appropriate 
then, 


w = T^{g{0 + T 2a -^W ! + ...), 

u _ r 7+*-0(_ Lyg' + T 1«-W Ux + ...), (58) 

P = Po(y,t) + T 2a ~ 2l3 P 1 (y,0+--; 

which on substitution into the axial momentum equation (26) and balance of terms yields, 

7 = /?-«, W ly = iy(P 04 - 9 "). 
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This in turn leads to the expression for U\ from the continuity equation and the result P y = 0. Use 
of these results into the tangential and normal stress balance equations (28)-(29) gives at 0(r~ f3 ) 
and 0(r~ a ) respectively, on y = /(£): 

(/V)' = \f 3 P0\ P 0 (t) = -g'+j, 

and on integration with respect to £ it yields 

1 k 

9 ~~Tf + J>' ( 59 ) 

with k a constant. A leading order balance in the kinematic condition (30) gives the second equation 

a = b (9 + POf' + (^9* ~ l)/= 0. (60) 

To leading order, then, the equations for the scaling functions (59) and (60) are identical to those 
found in Section 4.1 and solution proceeds as before, 

4.6 Validity of the similarity equations 

As the jet pinches, the following velocity field is found, 

W = T0- 1 g(Z) + T 1 -t } W 1 + ..., u= Ayg> + T 2-2/3 Ui+ ( 61 ) 

The leading order axial velocity becomes infinite as r ->• 0+ and since we are working within the 
context of the Stokes equations we need to check the effect of neglected inertial terms in the Navier- 
Stokes equations. More details may be found in Papageorgiou [28], [29]; briefly, it is found that the 
radial momentum equation approximation remains consistent as r — » 0; using the solutions (61) 
into the Navier-Stokes system shows that in deriving these solutions we keep terms of 0(t~ 1 ~ /3 ) 
while neglecting unsteady and nonlinear terms of 0(r^ 2 ). Consistency of the Stokes solutions all 
the way to pinching, is possible then if the neglected Navier-Stokes terms remain asymptotically 
smaller as r - 0. This requires fi - 2 < -1 - /?, or, fi > 1/2 (equality brings in unsteady and 
nonlinear terms - see below). It can be concluded, therefore, that the similarity solutions for Stokes 
jets which were found analytically and confirmed numerically, break down as the jet pinches and 
unsteady and inertial terms enter then. For high viscosity fluids one can expect the Stokes solutions 
to be operational for part of the evolution. It is particularly interesting to note that our analysis 
indicates that numerical algorithms based on the Stokes equations, as in the boundary integral 
method for example, need to be interpreted with care as the jet pinches. Our analysis quantifies 

exactly how such Stokes solutions may break down and in what follows we present the evolution 
that enters at that stage. 

5 Jet breakup for Navier-Stokes flows 

In this Section we will show how the breakup for Stokes jets needs to be modified in order to get a 
consistent solution. This is done in two ways: First we follow the analysis of Section 4.6. A system 
of similarity equations is found for the scaling functions. We also consider pinching according to 
the model equations of Section 3.1 with k ^ 0, and show that an identical system emerges. Finally 
we use these solutions to show how they can be used to continue the flow just beyond pinching by 
means of rational methods which obey the underlying physical principles, rather than introducing 
ad hoc regularizations as has been done by other authors. 
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5.1 Pinching solutions of the Navier-Stokes equations 

Guided by the estimates of the various terms as the Stokes jet pinches (see Section 4.6) we take 
/? = l/2 and look for terminal states of the Navier-Stokes equations of the following form: 


w = r“ 1/2 W 0 (0 + r l/2 Wi{y, f) + . . . , u = U 0 {y, f) + tU x { y, 0 + - • • , 

p = r~ 1 P(y,Z) + S(t,z) = r/(f), (62) 

r = ry, z = t 1/2 £, (63) 

where as before we take t = t„-t with t s the singular time and we assume, without loss of generality, 
that pinching occurs first at z = 0. The functions Wi, Ui, P are to be found. Substitution into the 
Navier-Stokes equations (25)-(26) gives to leading order, 

P y = 0, (64) 

Refywo + awfo + Wowt) = -P' + W Xyy + l -W Xy , (65) 

with primes denoting ^-derivatives. The solution for W x is 

W x = \y 2 ^ (P-W 0 + \Re(£W 0 + Wi)) + A(0 , (66) 


with A(£) some function of a solution for U\ follows from the continuity equation (27). A coupled 
system of equations for W 0 and / is found as follows: (i) evaluation of P from the leading order 
contributions to the normal stress balance (29) and insertion into (66) to give W\ y in terms of 
desired functions, (ii) substitution of Wi into the leading order contribution to the tangential stress 
balance equation (28) to yield one equation, (iii) evaluation of the leading terms of the kinematic 
equation (30) to give a second equation. The two equations are nonlinear ODEs and describe the 
s imil arity solutions within the pinch region. The equations read: 

3 Wg+£+ 6^ = R e Q(W> + eWS) + W)W') , (67) 

(Wo + \of'-^-\K)f = 0. ( 68 ) 

These equations are identical to those first derived and solved by Eggers [11], [12] by analysis of 
the model one- dimensional PDEs as briefly described next. 

5.2 Pinching solutions of the one-dimensional model equations 

The equations to be considered are the model system (40) and (41). Using the notation and 
methodology of previous Sections we look for pinching self-similar solutions in the form 

5 0 (M) = r/(0, «>o(<,s) = r-V 2 g(0, z = t 1 ' 2 £. (69) 

Substitution of (69) into (40) and (41) keeps all terms in the equations of the same order in r and 
yields exactly equations (67) and (68) where R e is now replaced by k; as mentioned earlier the 
constant k may be scaled out of the problem to arrive at exactly the system studied by Eggers [11], 
[ 12 ]. 

Numerical solutions of the initial value problem according to the model system (40) and (41) 
were first given by Eggers and Dupont [13] and also by Eggers [11]. It is found that as pinching is 
approached, general non-parity initial conditions evolve into universal scaling functions which are 
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not symmetric in the jet shape, for instance. Solution of the similarity equations (67) and (68) 
were obtained by Eggers [11], [12] by use of a shooting method. 

In our analysis of the dynamics beyond pinching, it is necessary to know the behavior of the 
scaling functions away from the pinch point, that is as the local inner solution is matched to the 
outer quasi-static solution. It is easily established that as |f| oo we have 

/(£W 2 , wb(o~r 1 . 

More precisely, there are universal constants bf and bf such that 

/(0~*f£ 2 , Wb(£) ~ bfr 1 as^^ioo. 

Numerical computation (see [12]) gives these constants to be 

b t = 4 - 6 35, b~ = 6.074 x 10~ 4 , 6+ = 0.0723, 6J = 57.043. 

Note that these forms at infinity imply that the jet shape is locally parabolic as a pinch forms 
unlike the corresponding shapes for Stokes jets. 


(70) 


5.3 Dynamics beyond pinching 

In this Section we follow the analysis of Papageorgiou [29] to show how the universal scaling 
functions described above can be used as the correct initial conditions for the flow just after the 
change of topology in order to obtain a theoretical description of the dynamics for small times just 

after pinching. The equivalent analysis for Stokes jets using the scaling functions of Section 4.2 can 
be found m [29]. 

r *5^ rt here ViS r C °? fl ° W analogUes of the methods originally developed by Taylor 

[43], Keller [18], Ting and Keller [44] and more recently Keller, King and Ting [19]. The solutions 
given in previous Sections describe the flow up to the time of breaking. Beyond this time the jet 
separates and the two ends accelerate away from the breaking point. Assuming that the jet ends at 
some axial position X{t) and that a fluid blob of mass M(t) is attached to that end, it is suggested 
that a fairly complete and rational picture of the dynamics beyond pinching is possible if one can 
solve for X(t) and M(t), for small times at least. These small time solutions (given below) are 
useful in initializing a direct numerical simulation correctly. We note that for the viscous case the 
exponents of the scaling functions are fixed and so the dynamics described below are universal. This 
is m contrast to the analysis of [44] where an exponent has to be chosen with different dynamics 
resulting. In addition this theory is rational in the sense that no ad hoc regularizations are made 
at can allow the jet to pinch and form fluid blobs even though the equations describing such 
motions are slender jet equations (see [12] for an ad hoc regularization). Instead we proceed by 

using momentum and mass balance arguments to obtain the evolution of X(t) and M(t) and use 
the correct initial conditions at t = 0+. 

The first equation comes from a mass balance of the fluid inside the blob. Define M(t) to be the 
blob mass and X(t) to be the axial distance from the breaking point where the blob is attached to 
the slender jet. The jet radius at the point of attachment is S(X(t),t) and W is the axial velocity 

there. All quantities are dimensional at this stage. The mass balance of fluid entering the blob 
through the attachment region is 


dM 

dt 



(71) 


Next consider balance of momentum using dimensional variables. The forces acting on the blob 
at * - X ( f ) are due to viscous stresses and the pull due to surface tension. Denote the axial 
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component of the viscous stress by T zz and the radial component by T rr , then the total force acting 
on the cross-section at X{t) is 

F = 7 rS 2 (T zz -p) + 2 izoS, 

which on use of the normal stress balance equation (for slender jets) 

_ rr , <r 
P — Trr "b 

becomes 

F = 7 tS 2 (T zz - T rr ) + X&S = Si zpS 2 W z + iraS. (72) 

Note that for Stokes jets F z is zero since Stokes flows do not possess inertia; this is recognized to 
be equation (40) with k = 0 (see [32] also). With inertia present, however, a balance of momentum 
of the fluid blob is found using (72): 

^ (m^P\ = SirpS 2 W z + izS + W^f. (73) 

dt \ dt ) dt 

On use of the capillary scales of Section 3 the dimensionless forms of equations (71) and (73) are 
found to be 

" = '* (£-"')• (74) 

r 4 ( m ?) = ( 75 > 

where the same symbols are used in denoting dimensionless variables. Note also that the inviscid 
system studied by Ting and Keller [44], has the viscous term SirS 2 W z term missing from (75). 

Equations (74) and (75) are to be solved subject to initial conditions, at t — 0 say, which are 
given by the terminal similarity states described in earlier Sections. In coupling the flow into the 
dynamics of the fluid blob, we require the behavior of the solutions for small but positive time. The 
Navier-Stokes were considered in Section 5 as the singular time was approached from below. The 
si mil arity equations just beyond pinching (assumed without loss of generality to occur at t = 0), 
in regions where the jet is slender at least, are obtained for 0 < t << 1 by the ansatz (62) with 
t > 0 replacing r; time derivatives are different now and following the same procedure the following 
equations need to be solved for the scaling functions / and g : 

f> w' f' 1 

3 Wg + + 6—y^- = -i(Wb + £W') + WoWi (76) 

(Wo-^)/ , + (l + ^o) = 0, (77) 

where R e = 1, since it can be scaled out of the problem. In Section 5 we described the universal 
solution as the jet breaks; in particular, then, equations (76) and (77) must be solved subject to the 
conditions at infinity given by (70). For instance, choosing to follow the dynamics of the rightmost 
fluid blob, we take the conditions for large positive £ provided by (70). These conditions are 
sufficient to determine unique scaling functions for smaller £ by integrating inwards from infinity. 
The computation can be valid only as far as the point where the fluid blob begins, i.e., at z = X(t) 
(the slender jet ansatz breaks down otherwise). Clearly then, we need to solve for X(t) in order to 
find the appropriate scaling functions which in turn provide the evolution of M(t). 
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Following Ting and Keller [44] we look for solutions of equations (74) and (75) for small t as 
follows: 


X(t) = X 0 t 1 ' 2 , S(t,X 0 ) = tf(X 0 ), W(t,X 0 ) = t- 1 / 2 W 0 (X 0 ). 

Substitution into (74)-(75) gives 


M(t) = ^tt/ 2 (X 0 ) Qx 0 - Wo(X 0 )) * 5/2 , (78) 

Qx 0 - W 0 (X 0 )j Qx 0 - W 0 (X 0 )j f(X o) - 3/(Xo)W'(Xo) = 1. (79) 

Note that equation. (79) provides a condition that determines Xo. Numerically this is done as 
follows: Starting at large positive £ and using the appropriate universal asymptotic forms (70) 
there, we integrate the system (76) and (77) inwards until the first point is found where equation 
(79) is satisfied. The numerical value found from the computations is Xo = 0.91 correct to two 
decimals. With Xo known (as well as the scaling functions there), we use (78) to compute the mass 
M(t) at different times. 

A schematic of the jet shapes, for instance, as time increases can be constructed as follows: We 
chose four different times t = 0.0001 — 0.0004. Xo is computed as described above and so are M(i) 
and X(t) Xo t ! . The point X(f) marks the end of the slender part of the jet and the beginning 

of the fluid blob. We picture the jet shape by distributing the mass M(t) over a sphere in order 
to compute its radius. Results of such calculations are given in [29]. We note that to within our 
treatment of the blob shape, the results are based on variables which are universal and so the small 
time solutions provide rational initial conditions for jet dynamics beyond breakup. 

6 Conclusions 

We have considered the breakup of liquid jets under the action of capillary instability. The phe- 
nomenon of drop formation which precedes the topological change in the shape has been analyzed 
by construction of local self-similar solutions for both viscous and inviscid jets. It has been demon- 
strated that viscous jets admit slender breakup and unique scaling exponents and functions. One- 
dimensional models based on a slender jet ansatz have also been solved numerically by specification 
of arbitrary smooth initial conditions, and exhibited to support the theoretical self-similar predic- 
tions. Leading order 1-D models are sufficient when viscosity is present, but higher order terms in 
the expansion must be incorporated into the evolution equations for inviscid flows. Such reduced 
models have also been analyzed for self-similar solutions as well as by numerical simulations, and 
the results indicate drop formation with a locally cusp-like behavior. The system governing the 
dynamics on the short scale that removes the cusp singularity is also given. 
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FIGURE 4 Stokes jets, 
parameter / 2 . 
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FIGURE 9 Simulation using 1-D Stokes model. Convergence to a scaling function as t — t s . 
(a) The scaling function /(£), (b) p(^). 
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FIGURE in Simulation using I-D Stokes model. Non-symmetric initial conditions. Evolution 
of S(z,t) and W{z,t). 
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FIGURE 11 Non-svmmetric initial conditions. 
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